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. Abstract. Our understanding of the Lyman-a forest has received a 

q_) ' great boost with the advent of the Keck Telescope and large 3D hydro- 

Q . dynamical simulations. We present new simulations using the SPH tech- 

nique with a P 3 MG (Particle-Particle Particle-MultiGrid) non-periodic 
gravity solver. Our method employs a high resolution (1 kpc) inner vol- 
ume, essential for capturing the complex gas physics, a larger low res- 
olution volume, essential for correct larger scale tidal fields and a self- 
consistently applied, uniform tidal field to model the influence of ultra 
long waves. We include a photoionizing UV flux and relevant atomic 
cooling processes. We use constrained field realisations to probe a selec- 
ts) ■ tion of environments and construct a statistical sample representative of 

the wider universe. We generate artificial Lyman-a spectra and fit Voigt 
■ profiles. We examine the importance of (1) the photoionizing flux level 

and history, (2) tidal environment and (3) differing cosmologies, including 



in 



oo 



■ CDM and CDM+A. With an appropriate choice for the UV flux, we find 



Oh 



that the data is fit quite well if the rms density contrast is ~ 1 at z ~ 3 



O ' on galaxy scales. 

■ To appear in Computational Astrophysics, Proceedings of the 12th Kingston Conference, 

Halifax, October 1996, ed. D. Clarke & M. West (PASP) 



1. Introduction 

We hope to be able to use the current wealth of Lyman-a absorption data to constrain 
the shape and normalisation of the density power spectrum and thus the cosmological 
model. The extraction of this information is complicated by gas physics, the ultraviolet 
flux, the star formation history and supernova energy injection. By employing gasdy- 
namical simulations, we can incorporate the detailed high-redshift environment of the 
clouds self-consistently. We perform our Lyman-a simulations using the SPH method 
in conjunction with a gravity solver based on the multigrid method. We describe our 
code, with attention to the novel features, in section ||] 

Simulations have to be carefully constructed to achieve both the necessary reso- 
lution and simultaneously provide a good sample of the universe. Large scale power is 
also of particular importance for producing the correct tidal env i ronment for Lyman-a 



clouds . Recent Lyman-a simulations (e.g. Rauch et al. (1996) 



Miralda-Escude et al. 



|(1996)| , iMiicket et al. (1996)| , [Dave et al. (1996)| and ^hang et al (1996)| ) have traded 



good fc-space resolution 210 — 17.95 h Mpc 1 kpc against size, 3 — 22 Mpc. We address 
these design issues in section m and see also Bond & Wadsley (1996), (hereafter BW). 
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In generating the line statistics of the artificial Lyman-a spect ra, it is important 
to fit the l i nes in a manner di rectly comparable with observations (Miralda-Escude ei 



aL (1996) , Dave et aL_(1996) ). Our method of fitting Voigt profiles is described with 



our results m section [4.| 



2. Numerical Method 



Smoothed Particle Hydrodynamics (SPH) is a fully Lagrangian method fo r fluid dy- 



namics . It has been demonstrated to be robust and flexible (e.g. Review by Monaghan 



SPH has the advantage of following collapse of structure with constant mass 
resolution. Eulerian codes have high gas resolution in voids, but structures arising 
from small scale perturbations in the early universe are limited to the mass resolution, 
which is similar in both methods. Eulerian codes without adaptive refinement smooth 
structure in collapsing regions. 

We use a predictor-corrector time stepping scheme allowing us to vary the global 
timestep. The code runs very quickly at first, but slows when the particle-particle 
section of the gravity solver becomes dominant. 

The code includes radiative cooling and photoionization heating with equilibrium 
abundances. The species we consider are H, H + , He, He + , He ++ and e~. 

We make use of a highly stable, 2nd order implicit scheme for the energy equation 
to avoid excessively small timesteps associated with the heating/cooling timescales. We 
find an iterative solution for the energy, E, 

E (t + At) - E{ t) pdv + _ / m+mt+AQ p 



At 

The predicted midpoint values for p and the PdV work are used. 
2.1. The P 3 MG Gravity Solver 

Iterative Full Multigrid techniques provide an excellent alternative to Fourier based 
PM Methods. Though comparable in speed for a periodic computation, the Multigrid 
method requires less work for free boundary calculations, by roughly a factor of two. In 
addition, we take advantage of previous time-step information to increase speed even 
further by reducing the number of iterations required for convergence. We employ a 
multipole expansion of the particle potentials to provide boundary conditions for the 
free boundary calculation. We have added Particle-Particle interactions to dramatically 
augment the resolution of our gravity solver. The resulting P 3 MG code achieves ~1% 
force errors. 

For these simulations we have multiple particle masses. We make use of nested 
grids to optimise our computational costs. We use a 128 3 grid with P-P around the 
high resolution core region with an additional 64 3 grid with no P-P forces for the lowest 
resolution particles whose function is only to provide long range tidal forces. This 
provides for a maximum resolution of 1 kpc. We set this as the resolution limit for our 
gas and gravity. 



2.2. Enhancements to Standard SPH 

SPH relics on weighted sums over near neighbours to estimate all fluid quantities in- 
cluding the pressure forces that move the particles. The method is adaptive, with the 
smoothing length for the weighted sum being a direct measure of the local resolution. 

Appropriately varying smoothing lengths, h, are essential to any SPH implementa- 
tion. Methods based on numerically estimated quantities such as divergence or density 
are susceptible to noise whereas what we require is a good sample of neighbouring parti- 
cles within 2 h. To this end we count neighbours in three equi- volume shells around each 
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Redshift Z Redshift Z 

Fi glirG 1. UV Flux history and Continuum Depression, Da- Curve (2) in the left 
hand panel was designed to reproduce the Da vs. redshift data shown in the right panel. 
A curve such as (3), including an estimate of the UV flux from early stars, is also consistent 
with the Da data. The flux model of Haardt & Madau (f996) is too high for standard 
CDM models, as demonstrated by the dashed curve in the right hand panel. 



particle, (0 — 1.58) h id, (1.58 — 2) h id and (2 — 2.29) h id- We update h by interpolating 
over these bins to have 36 neighbours within 2 h new . We perform the counting during 
the density summation (to avoid extra computation) for use during the next timestep. 

The SPH method involves only pairwise interactions, giving exact momentum and 
energy conservation. It also relies upon explicit artificial viscosity to model shocks 
correctly. The artificial viscosity for particle i, contributed by a neighbour j, is a 
function of (Uj — Vj ) • (fi — fj). This form cannot differentiate between truly converging 
flows and pure shear flows. Numerical viscosity is thus in discriminately appl ied unless 



a switch is used, such as \div\/(\div\ + \curl\ + e) (e.g. 3teinmetz (1996)). Rather 
than finding the exact divergence and vorticity to flag particles in strongly convergent 
(shocking) regions, we use a similar form that is easier to calculate, 

Ej Wijinj^i -vj)- (n - fj) 
Switch, = J — _ 7T7 • 

Lj w v m 3 I K -Vj)- (n -rj)\ 

A linkcd-list of particles allocated to a grid of cells is an extremely fast method for 
locating particle neighbours. Tree structure is an alternative method, though roughly 
six times slower for uniform densities and very memory intensive. However, with high 
density contrasts, many particles may accumulate in a few link list cells, rendering the 
method locally 0(n 2 ) and thus very slow. We adaptively refine our linked list grid 
to avoid this problem - to a maximum level of four binary sub-divisions. The storage 
requirements are naturally much larger but the resulting method slows by only a factor 
of three in extremely clustered arrangements. 



The Simulations 



The photoionization levels and heating are determined by the background UV flux, in 
particular the flux at the ionization edge for hydrogen, Jgi2- Observati ons of the proxim - 
ity effect hint at little evolution in Jgi2 over the range z = 2 to 4 (e.g. Bcchtold (1994)| ). 
As shown in the left panel of F igure fl we apply both a fixed flux (1) and a strongly 
evolving flux similar to that of Haardt fc Madau (1996) lowered to fit the continuum 
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depression data (2), given by the functional form, 5 x 10~ 22 ((l + z)/3) 9 exp(— 3 (z — 2)). 
These choices should roughly bracket the true mean UV flux. 

The small scale nature of the clouds demands high resolution numerical work, 
especially as the collisionless nature of the dark matter plus radiative cooling prevent 

the classical Jean's length, fcj 1 « 0.049(4/(1 + z) {T/10 4 K)) 1 / 2 n~r /2 h" 1 Mpc from 
setting a minimum scale. Q nr refers to density parameter in non-relativistic matter. 

We simulate regions of comoving diameter 5.0 Mpc (2.5h _1 Mpc for the CDM 
model, 3.5h _1 Mpc for the CDM+A model) with 50 3 /2 gas and 50 3 /2 dark matter 
particles at our highest resolution. We enclose this in a medium resolution buffer of 
gas and dark matter particles with 8 times the mass, out to 8.0 Mpc and then massive 
tide particles with 64 times the mass, out to 12.8 Mpc. We apply a linearly evolved 
external shear to the entire volume, derived from the initial conditions. Thus we make 
considerable effort to model the long wave, tidal environment well. 

We use a sophisticated constrained field code to set up the initial displacements 
from unperturbed lattice positions, using the Zeldovich approximation, ft allows for 
arbitrary types and numbers of constraints. Our initial conditions finely sample an 
enormous range of /c-space, down to k — 0.01 hMpc^ 1 for these simulations (see Figure 1 
of BW). Periodic simulations must truncate the large scale power beyond the scale of 
the fundamental mode of the box and the power spectrum is still very flat at k ~ 1 
hMpc -1 , the scale of Lyman-a simulations. 

Most simulations use random regions of the Universe. We prefer to have control 
parameters which govern some of the major large scale characteristics of the simulation 
volume. We constrain the density (via v), bulk flow and tidal fields smoothed over three 
Gaussian filtered scales Rf, where v — 5L, P k/cR f the overdensity of the region relative 
to the rms fluctuation level an f . More details and applications of our peak constraint 
approach are described in BW. For these simulations, we fix v for Rf — 0.5 Mpc and use 
mean field expectations to set the other parameters, and the values at two other scales, 
1.0 and 1.5 Mpc. The z/s used are shown in Table |l|. Rare events which form large 
"bright" galaxies by z ~ 3 in the patch require higher v. v has a Gaussian probability 
distribution, allowing us to appropriately weight and combine simulations of different v 
to create a sample broadly representative of the universe. In creating samples of Lyman- 
a statistics, the solid angle on the sky subtended and length of spectra generated by 
the region add further weighting. 



Table 1. Initial Conditions used in this work. 

Cosmology Description v e v , p v a Weight 6 Angle c A A c 

CDM Void XI 0.179, 0.014 O069 1747 OB" 

Weak Void -0.7 0.337, 0.034 0.242 1.26 1.36 

Flat 0.0 0, 0.379 1.0 1.0 

Weak Peak 0.7 0.337, 0.034 0.242 0.78 0.77 

Peak 1.4 0.179, 0.014 0.069 0.57 0.62 

CDM+A Flat 0.0 0, 



"Tidal (shear) field parameters 
6 Weighting due to v only 

c Solid angle and average length of spectra, A A, measured at z = 3 relative to Flat region 

We describe results for two cosmological models: standard CDM, normalised to 
er 8 = 0.67 with h = 0.5, Q B = 0.05; COBE-normalised CDM+A, with h = 0.7, tt B = 
0.0255 and VL nr = 0.335 (shape parameter T — 0.225 and tilt n s — 0.94, which gives 
a theory in agreement with the shape of the observed large scale power spectrum, see 
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Figure 1 of BW). We required that <7o.5Mpc = 1-05 at z = 3 and the number of baryons 
in the simulation volume to be the same in both models. 

4. Results 

To analyse our simulations we produce artificial spectra, with signal to noise of 100 and 
5 kms" 1 pixels. We sample a regular spatial grid of lines of sight running along the 
three axial directions through our simulations. A typical line of sight is ^500 kms -1 
long and a single simulation sample ~ 10 6 kms -1 at z=3. We fit Voigt profiles using 
an automated profile fitting program, designed to emulate the methods employed by 
observers. The program identifies each group of lines as a region in the spectrum where 
the flux drops below 98% of the continuum. \ 2 minimization is used to fit first one 
line, then two lines and so on. The number of lines that produces the lowest reduced- 
X 2 = X 2 /d is used, d is the number of degrees of freedom remaining for the fit, estimated 
as 1 per 2 pixels in the line group region minus 3 for each line used in the fit. In practice 
a line group rarely needs more than 6 lines for a good fit. 

Figure @ demonstrates visually the difference between several runs, as shown in the 
column depth of neutral hydrogen through the central 2 Mpc of the simulations. 

The left panels show the dramatic effect of varying v and hence the shear field. The 
filamentary structure is greatly enhanced in even slightly overdense regions. This has 
important repercussions in the types and numbers of Lya absorption lines produced, as 
shown in the bottom left panel of Figure ||. Simulating overdense regions or voids like 
these is not possible in a periodic box of a similar size, because the density averaged 
over the box must remain equal to the universal mean value. Even for significantly 
larger boxes this is a problem. 

Comparing the CDM and CDM+A runs (top left and right panels in Figure || 
respectively), it can be seen that the flatter spectrum in the CDM+A case makes the 
filaments more prominent and the dwarf galaxies less so, but without it having a major 
impact on the Njji frequency curve. The difference in normalisation is due to the 
width in rcdshift of the simulation volume in each cosmology, dz/dl com = iJo/c(O nr (l + 
z) 3 + ^a) 1 / 2 . The high resolution size is A/ com =5 Mpc, comoving, in all simulations. 
Dividing by the ratio of Az = dz/dl com A.l com at z=3 moves the CDM+A curve down 
by 0.09 on the top right panel of Figure [| 

The left hand panels of Figure demonstrate that an appropriate fixed choice of 
J912 = 5x 10~ 22 can reproduce the observed statistics very well. We have statistically 
combined 5 CDM simulations to effectively produce a good sample of the universe that 
includes rarer regions. The advantage of this over a single large simulation is great 
resolution and importance sampling. The contribution of different regions is apparent. 
The voids, v < 0, cause a lowering of the curve, especially at low column depths. 

The results of a low resolution run (half that of the standard ones) are shown 
represented by crosses in the top right panel of Figure ||[ Though it underpredicts 
the low column lines, the high column lines are overpredicted relative to our higher 
resolution runs: limited resolution mimicks heating by "puffing" up the cores of cold 
lumps out to greater sizes and thus gives higher cross sections. 

In the top right panel of Figure || we demonstrate that attempts to renormalise 
simulations with inappropriate values of the ionizing flux will fail to give the correct 
frequency of Njji- The curves of interest are the CDM line with a fixed UV history 
and the scaled low UV points (triangles). These simulations are otherwise identical and 
have been processed with a fixed value of J912 = 5 x 10 -22 . What does scale fairly 
well is the continuum depression, Da, shown in the right panel of Figure [|. Da is 
defined at the mean value of 1 — exp(— r) for a spectral region. This is because the 
dominant contributors to the optical depth are lines with Njji ~ 10 14 , resulting from 
barely overdense material for which the UV flux history is less important. 
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Figure 2. HI Column depth through a cube 2 comoving Mpc thick at z=3. The 
contour levels, in logician" 2 ), are 11,11.25,11.5,11.75,12,12.25,12.5,12.75,13 (dotted) ,13.36 
(mean density, dot dashed) and 14,14.5,15,17 (solid). The left panels are CDM simulations, 
with v = 0, +1.4 and —1.4, top to bottom. Note the incredible enhancement of the 
filamentary structure in the overdense regions. The right hand panels are v = runs 
directly comparable with the CDM run in the top left panel. In the top right panel we 
show a CDM+A run, which has more large scale power and thus more enhanced filaments. 
The run in the right centre panel had a lower UV flux at early times, resulting in a more 
clumpy distribution and greatly enhanced HI column depths. The lower right panel shows 
that poor resolution (1/2 the spatial resolution of the upper left panel) will remove much 
of the structure. 
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Figure 3. Frequency distribution of lines by Njji- In all frames the boxes represent 
binned data, box width equalling bin width and the upper and lower edges denoting 1 
sigma Poisson errors. At the top left the full range in frequency is plotted and in the 
remaining graphs, the power law — 1.46 logio(Nfij) + 8.25 has been subtracted for clarity. 
On the left, the data from 5 CDM simulations are plotted. The symbols in the lower plot 
represent (from the top down) the v = 1.4, 0.7, 0, —0.7 and —1.4 runs respectively and the 
combined data are shown with 1 sigma Poisson error bars. At the top right, the effects of 
various choices for the UV flux, numerical resolution and cosmology are compared. The 
low UV flux at early times results in more dumpiness in the IGM, overproducing lines at 
intermediate column depths. At the bottom right, the data resulting from combining the 
CDM simulations arc shown at 3 rcdshifts, divided by ((1 + z)/4) 2 . The evolution inferred 
from the observations is closer to (1 + z) 3 , suggesting a UV flux decreasing with redshift 
may be needed. 
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The past history of the UV flux is important for the formation of dense gaseous 
clumps, as illustrated in the top left and centre right panels of Figure @; the same 
UV flux was applied for the post-processing analysis, but the simulation shown in the 
centre right panel experienced a low UV flux during its evolution (curve (2) in Figure [l]). 
Applying a different flux in post-processing the calculation can of course have no effect 
on the density distribution and thus cannot co mpensate for a d ifferent history. 



Subtracting a power law, derived from the Hu et al. (1995) low end slope, from the 
f{Njji) curve of Figure reveals interesting features in the data. If the bumps at Nui — 
10 14 and N HI = 10 17 (where self shielding of the UV flux starts to become important) 
are real, they may allow simulations to finely discriminate between parameter choices. 
The data appear to follow a power law slope of -1.9 in the range Nhi = 10 14 — 10 16 cm~ 2 , 
which our runs with fixed flux reproduce and t hose with low UV flux at early times do 



not. This is also apparent in the results of |Rauch et al. (1996) who apply a UV 



flux that slowly increases with rcdshift. Simulation s that calculate the UV flux history 



(Miralda-Escude et al. (1996), Miicket et al. (1996)) perform fairly well, suggesting that 
a more detailed self-consistent treatment could produce a strong result. The excess of 
intermediate column lines seen for the low UV simulations may be due to unrealistically 
low energy input at early times; specifically there should be input from early stars. 
We are testing this hypothesis using a UV flux history similar to curve (3) in Figure [j]. 
Star formation and supernova feedback would also inhibit the formation of the extremely 
dense, cool clumps associated with this excess. We are now incorporating these processes 
using a general chemical evolution code. 

We have illustrated that the Lyman alpha forest is a sensitive probe of the conditions 
that operated at z ~ 5 to 2. It is clear that we must produce a comprehensive sample 
if we are to deliver reliable predictions to the observers of the clouds. Finding forest 
observables that allow us to separate the influence of primordial power spectrum shape 
and amplitude, cosmological parameters such as fl^, Q c dm, Q-hdmi Hq and Qb, the 
ionization history and local energy injection from galaxies (feedback) will be a challenge 
to the growing group of cloud simulators, in close collaboration with the observers of 
the forest. Nonetheless it is clear that the basic model ( whether CDM or CDM+K ) with 
Co.5Mpc ~ 1 at z ~ 3 provides a rather good fit to the data. 
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